clc;
clear;

%% main主脚本
tic;
startTime = cputime;
pde = process_pde_HMF_t(1);
vector_h = [1/2 1/4 1/8 1/16];
basis_type1 = "P2";
basis_type2 = "P1";

u_Linf = zeros(size(vector_h,2),1);
u_L2 = zeros(size(vector_h,2),1);
u_H1 = zeros(size(vector_h,2),1);
w_Linf = zeros(size(vector_h,2),1);
w_L2 = zeros(size(vector_h,2),1);
w_H1 = zeros(size(vector_h,2),1);
T_Linf = zeros(size(vector_h,2),1);
T_L2 = zeros(size(vector_h,2),1);
T_H1 = zeros(size(vector_h,2),1);
p_Linf = zeros(size(vector_h,2),1);
p_L2 = zeros(size(vector_h,2),1);
p_H1 = zeros(size(vector_h,2),1);

for i = 1:size(vector_h,2)
    h = vector_h(i);
    dt = h^2;
    Nt = ceil((pde.domain(6) - pde.domain(5))/dt);

    mesh = process_mesh(pde, h);
    fs1 = process_fespace(mesh, basis_type1);
    fs2 = process_fespace(mesh, basis_type2);

    % process_showmesh(mesh, fs1);
    % process_showmesh(mesh, fs2);

    A11 = term_U_V(mesh, fs1);
    A13 = (pde.nu+pde.nu_r)*term_grad_U_grad_V(mesh, fs1);
    A14 = term_P_div_V(mesh, fs1, fs2);
    A15 = 2*pde.nu_r*term_rot_W_V(mesh, fs1);
    A16 = pde.eh*term_JT_V(mesh, fs1);

    A31 = term_W_Psi(mesh, fs1);
    A33 = (pde.Ca+pde.Cd)*term_grad_W_grad_Psi(mesh, fs1);
    A34 = 4*pde.nu_r*A31;
    A35 = 2*pde.nu_r*term_rot_U_Psi(mesh, fs1);

    A41 = A31;
    A43 = pde.ka*term_grad_T_grad_G(mesh, fs1);

    O1 = sparse(fs1.Nb, fs1.Nb);
    O2 = sparse(fs1.Nb, fs2.Nb);
    O3 = sparse(fs1.Nb, 2*fs1.Nb);
    O4 = sparse(fs2.Nb, fs2.Nb);
    O5 = sparse(fs2.Nb, 1);

    A = [1/dt*A11+A13, -A15, -A16, -A14;
        -A35, 1/dt*A31+A33+A34, O1, O2;
        O3, O1, 1/dt*A41+A43, O2;
        -A14', O2', O2',O4];

    U1 = [pde.u1(fs1.node(:,1), fs1.node(:,2), 0);
        pde.u2(fs1.node(:,1), fs1.node(:,2), 0)];
    W1 = pde.w3(fs1.node(:,1), fs1.node(:,2), 0);
    T1 = pde.T(fs1.node(:,1), fs1.node(:,2), 0);

    for m = 0:Nt-1
        t = (m+1)*dt;
        A12 = term_Uold_grad_U_V(mesh, fs1, U1);
        A32 = term_Uold_grad_W_Psi(mesh, fs1, U1);
        A42 = term_Uold_grad_T_G(mesh, fs1, U1);
        A44 = pde.D*term_rot_Wold_grad_T_G(mesh, fs1, W1);
        AN = blkdiag(A12, A32, A42-A44, O4);
        At = A + AN;

        b1 = term_F_V_t(mesh, fs1, pde.f11, pde.f12, t);
        b2 = term_F_Psi_t(mesh, fs1, pde.f2, t);
        b3 = term_F_Psi_t(mesh, fs1, pde.f3, t);
        bt = [b1 + 1/dt*A11*U1;
            b2 + 1/dt*A31*W1;
            b3 + 1/dt*A41*T1;
            O5];

        [At, bt] = process_boundary_HMF_t(At, bt, mesh, fs1, fs2, pde, t);

        SOL = At\bt;

        U1 = SOL(1:2*fs1.Nb);
        W1 = SOL(2*fs1.Nb+1:3*fs1.Nb);
        T1 = SOL(3*fs1.Nb+1:4*fs1.Nb);
    end

    % process_plot_scalar(mesh, fs2, SOL(3*fs1.Nb+1:end));

    [u_Linf(i), w_Linf(i), T_Linf(i), p_Linf(i)] = process_error_HMF_t(mesh, fs1, fs2, pde, SOL, "Linf");
    [u_L2(i), w_L2(i), T_L2(i), p_L2(i)] = process_error_HMF_t(mesh, fs1, fs2, pde, SOL, "L2");
    [u_H1(i), w_H1(i), T_H1(i), p_H1(i)] = process_error_HMF_t(mesh, fs1, fs2, pde, SOL, "H1");

    fprintf("已完成\t%d/%d\n", i, size(vector_h,2));
end

process_rate(u_Linf, vector_h);
process_rate(u_L2, vector_h);
process_rate(u_H1, vector_h);
process_rate(w_Linf, vector_h);
process_rate(w_L2, vector_h);
process_rate(w_H1, vector_h);
process_rate(T_Linf, vector_h);
process_rate(T_L2, vector_h);
process_rate(T_H1, vector_h);
process_rate(p_Linf, vector_h);
process_rate(p_L2, vector_h);
process_rate(p_H1, vector_h);

TOCtime = toc;
CPUtime = cputime - startTime;
disp(["TOC time(s)", num2str(TOCtime)]);
disp(["CPU time(s)", num2str(CPUtime)]);



%% process_pde_HMF_t
function pde = process_pde_HMF_t(index)
switch index
    case 1
        pde.index = index;
        nu = 1;
        nu_r = 1;
        eh = 1;
        C0 = 1;
        Ca = 1;
        Cd = 1;
        ka = 1;
        D = 1;
        pde.nu = nu;
        pde.nu_r = nu_r;
        pde.eh = eh;
        pde.C0 = C0;
        pde.Ca = Ca;
        pde.Cd = Cd;
        pde.ka = ka;
        pde.D = D;

        pde.u1 = @(x,y,t)pi*sin(pi*x).^2.*sin(2*pi*y).*sin(t);
        pde.u2 = @(x,y,t)-pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        pde.w3 = @(x,y,t)pi*sin(pi*x).^2.*sin(2*pi*y).*sin(t)-pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        pde.T = @(x,y,t)pi*sin(pi*x).^2.*sin(2*pi*y).*sin(t)+pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        pde.p = @(x,y,t)10*cos(pi*x).*sin(pi*y).*sin(t);

        pde.u1dx = @(x,y,t)pi.^2.*cos(x.*pi).*sin(x.*pi).*sin(y.*pi.*2.0).*sin(t).*2.0;
        pde.u1dy = @(x,y,t)pi.^2.*cos(y.*pi.*2.0).*sin(x.*pi).^2.*sin(t).*2.0;

        pde.u2dx = @(x,y,t)pi.^2.*cos(x.*pi.*2.0).*sin(y.*pi).^2.*sin(t).*-2.0;
        pde.u2dy = @(x,y,t)pi.^2.*cos(y.*pi).*sin(x.*pi.*2.0).*sin(y.*pi).*sin(t).*-2.0;

        pde.w3dx = @(x,y,t)pi.^2.*sin(pi.*(x.*2.0-y)).*sin(y.*pi).*sin(t).*2.0;
        pde.w3dy = @(x,y,t)pi.^2.*sin(x.*pi).*sin(pi.*(x-y.*2.0)).*sin(t).*2.0;

        pde.Tdx = @(x,y,t)pi.^2.*sin(y.*pi).*sin(pi.*(x.*2.0+y)).*sin(t).*2.0;
        pde.Tdy = @(x,y,t)pi.^2.*sin(x.*pi).*sin(pi.*(x+y.*2.0)).*sin(t).*2.0;

        pde.pdx = @(x,y,t)pi.*sin(x.*pi).*sin(y.*pi).*sin(t).*-1.0e+1;
        pde.pdy = @(x,y,t)pi.*cos(x.*pi).*cos(y.*pi).*sin(t).*1.0e+1;

        pde.f11 = @(x,y,t)pi.*sin(x.*pi).^2.*sin(y.*pi.*2.0).*cos(t)-pi.*sin(x.*pi).*sin(y.*pi).*sin(t).*1.0e+1-pi.^3.*sin(y.*pi.*2.0).*sin(t).*(nu+nu_r).*(cos(x.*pi.*2.0).*2.0-1.0).*2.0-nu_r.*pi.^2.*sin(x.*pi).*sin(pi.*(x-y.*2.0)).*sin(t).*4.0+pi.^3.*cos(x.*pi).*sin(x.*pi).^3.*sin(y.*pi.*2.0).^2.*sin(t).^2.*2.0-pi.^3.*cos(y.*pi.*2.0).*sin(x.*pi).^2.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*sin(t).^2.*2.0;
        pde.f12 = @(x,y,t)-pi.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*cos(t)-eh.*pi.*sin(t).*(sin(x.*pi.*2.0).*sin(y.*pi).^2+sin(x.*pi).^2.*sin(y.*pi.*2.0))+pi.*cos(x.*pi).*cos(y.*pi).*sin(t).*1.0e+1+pi.^3.*sin(x.*pi.*2.0).*sin(t).*(nu+nu_r).*(cos(y.*pi.*2.0).*2.0-1.0).*2.0+nu_r.*pi.^2.*sin(pi.*(x.*2.0-y)).*sin(y.*pi).*sin(t).*4.0+pi.^3.*cos(y.*pi).*sin(x.*pi.*2.0).^2.*sin(y.*pi).^3.*sin(t).^2.*2.0-pi.^3.*cos(x.*pi.*2.0).*sin(x.*pi).^2.*sin(y.*pi).^2.*sin(y.*pi.*2.0).*sin(t).^2.*2.0;
        pde.f2 = @(x,y,t)pi.^3.*sin(t).*(Ca+Cd).*(sin(x.*pi.*2.0).*(-1.0./2.0)+sin(y.*pi.*2.0)./2.0+sin(pi.*(x-y).*2.0)).*4.0+nu_r.*pi.^2.*sin(t).*(sin(x.*pi).^2.*sin(y.*pi).^2.*-4.0+sin(x.*pi).^2+sin(y.*pi).^2).*4.0-pi.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*cos(t)+pi.*sin(x.*pi).^2.*sin(y.*pi.*2.0).*cos(t)-nu_r.*pi.*sin(t).*(sin(x.*pi.*2.0).*sin(y.*pi).^2-sin(x.*pi).^2.*sin(y.*pi.*2.0)).*4.0-pi.^3.*sin(x.*pi).*sin(x.*pi.*2.0).*sin(y.*pi).^2.*sin(pi.*(x-y.*2.0)).*sin(t).^2.*2.0+pi.^3.*sin(pi.*(x.*2.0-y)).*sin(x.*pi).^2.*sin(y.*pi).*sin(y.*pi.*2.0).*sin(t).^2.*2.0;
        pde.f3 = @(x,y,t)pi.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*cos(t)+pi.*sin(x.*pi).^2.*sin(y.*pi.*2.0).*cos(t)+pi.^3.*cos(x.*pi).*sin(x.*pi).^3.*sin(y.*pi.*2.0).^2.*sin(t).^2.*2.0-pi.^3.*cos(y.*pi).*sin(x.*pi.*2.0).^2.*sin(y.*pi).^3.*sin(t).^2.*2.0-ka.*pi.^3.*cos(x.*pi).^2.*sin(y.*pi.*2.0).*sin(t).*2.0-ka.*pi.^3.*cos(y.*pi).^2.*sin(x.*pi.*2.0).*sin(t).*2.0+ka.*pi.^3.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*sin(t).*6.0+ka.*pi.^3.*sin(x.*pi).^2.*sin(y.*pi.*2.0).*sin(t).*6.0+pi.^3.*cos(x.*pi.*2.0).*sin(x.*pi).^2.*sin(y.*pi).^2.*sin(y.*pi.*2.0).*sin(t).^2.*2.0-pi.^3.*cos(y.*pi.*2.0).*sin(x.*pi).^2.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*sin(t).^2.*2.0+D.*pi.^4.*cos(x.*pi).^2.*cos(y.*pi).^2.*sin(x.*pi).^2.*sin(y.*pi).^2.*sin(t).^2.*3.2e+1-D.*pi.^4.*cos(x.*pi.*2.0).*cos(y.*pi.*2.0).*sin(x.*pi).^2.*sin(y.*pi).^2.*sin(t).^2.*8.0;
        pde.domain = [0 1 0 1 0 1];
    case 2
        pde.index = index;
        nu = 1;
        nu_r = 1;
        eh = 1;
        C0 = 1;
        Ca = 1;
        Cd = 1;
        ka = 1;
        D = 1;
        pde.nu = nu;
        pde.nu_r = nu_r;
        pde.eh = eh;
        pde.C0 = C0;
        pde.Ca = Ca;
        pde.Cd = Cd;
        pde.ka = ka;
        pde.D = D;

        pde.u1 = @(x,y,t)10*x.^2.*(x-1).^2.*y.*(y-1).*(2*y-1).*cos(t);
        pde.u2 = @(x,y,t)-10*x.*(x-1).*(2*x-1).*y.^2.*(y-1).^2.*cos(t);
        pde.w3 = @(x,y,t)10*x.^2.*(x-1).^2.*y.*(y-1).*(2*y-1).*cos(t)+10*x.*(x-1).*(2*x-1).*y.^2.*(y-1).^2.*cos(t);
        pde.T = @(x,y,t)10*x.^2.*(x-1).^2.*y.*(y-1).*(2*y-1).*cos(t)-10*x.*(x-1).*(2*x-1).*y.^2.*(y-1).^2.*cos(t);
        pde.p = @(x,y,t)10*(2*x-1).*(2*y-1).*cos(t);

        pde.u1dx = @(x,y,t)x.*y.*cos(t).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+1;
        pde.u1dy = @(x,y,t)x.^2.*cos(t).*(x-1.0).^2.*(y.*-6.0+y.^2.*6.0+1.0).*1.0e+1;

        pde.u2dx = @(x,y,t)y.^2.*cos(t).*(y-1.0).^2.*(x.*-6.0+x.^2.*6.0+1.0).*-1.0e+1;
        pde.u2dy = @(x,y,t)x.*y.*cos(t).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*-2.0e+1;

        pde.w3dx = @(x,y,t)y.*cos(t).*(y-1.0).*(x.*2.0+y-x.^2.*y.^2.*6.0-x.*y.*1.0e+1+x.*y.^2.*6.0+x.^2.*y.*1.8e+1-x.^3.*y.*8.0-x.^2.*6.0+x.^3.*4.0-y.^2).*-1.0e+1;
        pde.w3dy =  @(x,y,t)x.*cos(t).*(x-1.0).*(x+y.*2.0-x.^2.*y.^2.*6.0-x.*y.*1.0e+1+x.*y.^2.*1.8e+1+x.^2.*y.*6.0-x.*y.^3.*8.0-x.^2-y.^2.*6.0+y.^3.*4.0).*-1.0e+1;

        pde.Tdx = @(x,y,t)y.*cos(t).*(y-1.0).*(x.*2.0-y+x.^2.*y.^2.*6.0+x.*y.*2.0-x.*y.^2.*6.0+x.^2.*y.*6.0-x.^3.*y.*8.0-x.^2.*6.0+x.^3.*4.0+y.^2).*-1.0e+1;
        pde.Tdy = @(x,y,t)x.*cos(t).*(x-1.0).*(-x+y.*2.0+x.^2.*y.^2.*6.0+x.*y.*2.0+x.*y.^2.*6.0-x.^2.*y.*6.0-x.*y.^3.*8.0+x.^2-y.^2.*6.0+y.^3.*4.0).*1.0e+1;

        pde.pdx = @(x,y,t)cos(t).*(y.*2.0-1.0).*2.0e+1;
        pde.pdy = @(x,y,t)cos(t).*(x.*2.0e+1-1.0e+1).*2.0;

        pde.f11 = @(x,y,t)cos(t).*(y.*4.0e+1-2.0e+1)-cos(t).*(y.*2.0-1.0).*(nu+nu_r).*(-y+x.^2.*y.^2.*6.0+x.*y.*6.0-x.*y.^2.*6.0-x.^2.*y.*6.0+x.^2.*3.0-x.^3.*6.0+x.^4.*3.0+y.^2).*2.0e+1+nu_r.*x.*cos(t).*(x-1.0).*(x+y.*2.0-x.^2.*y.^2.*6.0-x.*y.*1.0e+1+x.*y.^2.*1.8e+1+x.^2.*y.*6.0-x.*y.^3.*8.0-x.^2-y.^2.*6.0+y.^3.*4.0).*2.0e+1-x.^2.*y.*sin(t).*(y.*2.0-1.0).*(x-1.0).^2.*(y-1.0).*1.0e+1-x.^3.*y.^2.*cos(t).^2.*(x.*2.0-1.0).*(x-1.0).^3.*(y-1.0).^2.*(y.*-6.0+y.^2.*6.0+1.0).*1.0e+2+x.^3.*y.^2.*cos(t).^2.*(y.*2.0-1.0).*(x-1.0).^2.*(y-1.0).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+2;
        pde.f12 = @(x,y,t)cos(t).*(x.*4.0e+1-2.0e+1)+cos(t).*(x.*2.0-1.0).*(nu+nu_r).*(-x+x.^2.*y.^2.*6.0+x.*y.*6.0-x.*y.^2.*6.0-x.^2.*y.*6.0+x.^2+y.^2.*3.0-y.^3.*6.0+y.^4.*3.0).*2.0e+1-nu_r.*y.*cos(t).*(y-1.0).*(x.*2.0+y-x.^2.*y.^2.*6.0-x.*y.*1.0e+1+x.*y.^2.*6.0+x.^2.*y.*1.8e+1-x.^3.*y.*8.0-x.^2.*6.0+x.^3.*4.0-y.^2).*2.0e+1+x.*y.^2.*sin(t).*(x.*2.0-1.0).*(x-1.0).*(y-1.0).^2.*1.0e+1-x.^2.*y.^3.*cos(t).^2.*(y.*2.0-1.0).*(x-1.0).^2.*(y-1.0).^3.*(x.*-6.0+x.^2.*6.0+1.0).*1.0e+2-eh.*x.*y.*cos(t).*(x-1.0).*(y-1.0).*(x-y-x.*y.^2.*2.0+x.^2.*y.*2.0-x.^2+y.^2).*1.0e+1+x.^2.*y.^3.*cos(t).^2.*(x.*2.0-1.0).*(x-1.0).*(y-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+2;
        pde.f2 = @(x,y,t)-cos(t).*(Ca+Cd).*(x.*2.0e+1+y.*2.0e+1-x.^2.*y.^2.*7.2e+2+x.^2.*y.^3.*2.4e+2+x.^3.*y.^2.*2.4e+2-x.*y.*2.4e+2+x.*y.^2.*6.0e+2+x.^2.*y.*6.0e+2-x.*y.^3.*4.8e+2-x.^3.*y.*4.8e+2+x.*y.^4.*1.2e+2+x.^4.*y.*1.2e+2-x.^2.*1.2e+2+x.^3.*1.6e+2-x.^4.*6.0e+1-y.^2.*1.2e+2+y.^3.*1.6e+2-y.^4.*6.0e+1)+nu_r.*cos(t).*(x.^2.*y.^2.*1.2e+2-x.^2.*y.^3.*1.2e+2-x.^3.*y.^2.*1.2e+2+x.^2.*y.^4.*6.0e+1+x.^4.*y.^2.*6.0e+1-x.*y.^2.*6.0e+1-x.^2.*y.*6.0e+1+x.*y.^3.*1.2e+2+x.^3.*y.*1.2e+2-x.*y.^4.*6.0e+1-x.^4.*y.*6.0e+1+x.^2.*1.0e+1-x.^3.*2.0e+1+x.^4.*1.0e+1+y.^2.*1.0e+1-y.^3.*2.0e+1+y.^4.*1.0e+1).*2.0-x.*y.^2.*sin(t).*(x.*2.0-1.0).*(x-1.0).*(y-1.0).^2.*1.0e+1-x.^2.*y.*sin(t).*(y.*2.0-1.0).*(x-1.0).^2.*(y-1.0).*1.0e+1+x.^2.*y.^2.*cos(t).^2.*(x.*2.0-1.0).*(x-1.0).^2.*(y-1.0).^2.*(x+y.*2.0-x.^2.*y.^2.*6.0-x.*y.*1.0e+1+x.*y.^2.*1.8e+1+x.^2.*y.*6.0-x.*y.^3.*8.0-x.^2-y.^2.*6.0+y.^3.*4.0).*1.0e+2-x.^2.*y.^2.*cos(t).^2.*(y.*2.0-1.0).*(x-1.0).^2.*(y-1.0).^2.*(x.*2.0+y-x.^2.*y.^2.*6.0-x.*y.*1.0e+1+x.*y.^2.*6.0+x.^2.*y.*1.8e+1-x.^3.*y.*8.0-x.^2.*6.0+x.^3.*4.0-y.^2).*1.0e+2+nu_r.*x.*y.*cos(t).*(x-1.0).*(y-1.0).*(x+y-x.*y.*4.0+x.*y.^2.*2.0+x.^2.*y.*2.0-x.^2-y.^2).*4.0e+1;
        pde.f3 = @(x,y,t)ka.*cos(t).*(x-y-x.^2.*y.^3.*1.2e+1+x.^3.*y.^2.*1.2e+1-x.*y.^2.*6.0+x.^2.*y.*6.0+x.*y.^4.*6.0-x.^4.*y.*6.0-x.^3.*4.0+x.^4.*3.0+y.^3.*4.0-y.^4.*3.0).*2.0e+1+x.*y.^2.*sin(t).*(x.*2.0-1.0).*(x-1.0).*(y-1.0).^2.*1.0e+1-x.^2.*y.*sin(t).*(y.*2.0-1.0).*(x-1.0).^2.*(y-1.0).*1.0e+1-D.*x.^2.*y.^2.*cos(t).^2.*(x-1.0).^2.*(y-1.0).^2.*(x.*-1.0e+1-y.*1.0e+1+x.^2.*y.^2.*2.8e+1+x.*y.*2.8e+1-x.*y.^2.*2.8e+1-x.^2.*y.*2.8e+1+x.^2.*1.0e+1+y.^2.*1.0e+1+3.0).*2.0e+2-x.^2.*y.^2.*cos(t).^2.*(x.*2.0-1.0).*(x-1.0).^2.*(y-1.0).^2.*(-x+y.*2.0+x.^2.*y.^2.*6.0+x.*y.*2.0+x.*y.^2.*6.0-x.^2.*y.*6.0-x.*y.^3.*8.0+x.^2-y.^2.*6.0+y.^3.*4.0).*1.0e+2-x.^2.*y.^2.*cos(t).^2.*(y.*2.0-1.0).*(x-1.0).^2.*(y-1.0).^2.*(x.*2.0-y+x.^2.*y.^2.*6.0+x.*y.*2.0-x.*y.^2.*6.0+x.^2.*y.*6.0-x.^3.*y.*8.0-x.^2.*6.0+x.^3.*4.0+y.^2).*1.0e+2;
        pde.domain = [0 1 0 1 0 1];
    otherwise
        error("Invalid pde index.");
end
fprintf("pde index: %d\n", pde.index);
fprintf("pde domain: xmin,xmax,ymin,ymax\n");
disp(pde.domain);
end

%% term_JT_V
function matrix_A = term_JT_V(mesh, fespace)
Nb = fespace.Nb;
A1 = sparse(Nb, Nb);
A2 = term_P_Q(mesh, fespace);
matrix_A = [A1; A2];
end

%% term_grad_T_grad_G
function matrix_A = term_grad_T_grad_G(mesh, fespace)
matrix_A = term_grad_P_grad_Q(mesh, fespace);
end

%% term_W_Psi
function matrix_A = term_W_Psi(mesh, fespace)
matrix_A = term_P_Q(mesh, fespace);
end

%% term_gradW_gradPsi
function matrix_A = term_grad_W_grad_Psi(mesh, fespace)
matrix_A = term_grad_P_grad_Q(mesh, fespace);
end

%% term_Uold_grad_T_G
function matrix_A = term_Uold_grad_T_G(mesh, fespace, Uold)
matrix_A = term_Uold_grad_W_Psi(mesh, fespace, Uold);
end

%% rocess_boundary_HMF_t
function [A, b] = process_boundary_HMF_t(A, b, mesh, fespace1, fespace2, pde, t)
Nb1 = fespace1.Nb;
% Dirichlet
for k = 1:4
    % u
    bdn1 = bdary_node(mesh, fespace1, k);
    A(bdn1,:) = 0;
    A(sub2ind(size(A), bdn1, bdn1)) = 1;
    b(bdn1) = pde.u1(fespace1.node(bdn1,1),fespace1.node(bdn1,2),t);
    bdn2 = bdn1 + Nb1;
    A(bdn2,:) = 0;
    A(sub2ind(size(A), bdn2, bdn2)) = 1;
    b(bdn2) = pde.u2(fespace1.node(bdn1,1),fespace1.node(bdn1,2),t);
    % w
    bdn2 = bdn1 + 2*Nb1;
    A(bdn2,:) = 0;
    A(sub2ind(size(A), bdn2, bdn2)) = 1;
    b(bdn2) = pde.w3(fespace1.node(bdn1,1),fespace1.node(bdn1,2),t);
    % T
    bdn2 = bdn1 + 3*Nb1;
    A(bdn2,:) = 0;
    A(sub2ind(size(A), bdn2, bdn2)) = 1;
    b(bdn2) = pde.T(fespace1.node(bdn1,1),fespace1.node(bdn1,2),t);
end
% stablization
M = term_P_Q(mesh, fespace2);
matrix_A = blkdiag(sparse(4*Nb1,4*Nb1), -1e-6*M);
A = A + matrix_A;
end

%% process_error_HMF_t
function [error_u, error_w, error_T, error_p] = process_error_HMF_t(mesh, fespace1, fespace2, pde, SOL, error_type)
Nb1 = fespace1.Nb;
t_end = pde.domain(6);
switch error_type
    case "Linf"
        err_u1 = 0;
        err_u2 = 0;
        err_w3 = 0;
        err_T = 0;
        err_p = 0;
        for n = 1:mesh.Ne
            V = mesh.node(mesh.elem(n,:),:);
            [int_points_loc, ~] = integral_loc(V, "2d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "x");
            psi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "x");
            u1h = phi * SOL(fespace1.elem(n,:));
            u2h = phi * SOL(fespace1.elem(n,:)+Nb1);
            w3h = phi * SOL(fespace1.elem(n,:)+2*Nb1);
            Th = phi * SOL(fespace1.elem(n,:)+3*Nb1);
            ph = psi * SOL(fespace2.elem(n,:)+4*Nb1);
            f1 = max(abs(pde.u1(int_points_loc(:,1),int_points_loc(:,2),t_end) - u1h));
            f2 = max(abs(pde.u2(int_points_loc(:,1),int_points_loc(:,2),t_end) - u2h));
            f3 = max(abs(pde.w3(int_points_loc(:,1),int_points_loc(:,2),t_end) - w3h));
            f4 = max(abs(pde.T(int_points_loc(:,1),int_points_loc(:,2),t_end) - Th));
            f5 = max(abs(pde.p(int_points_loc(:,1),int_points_loc(:,2),t_end) - ph));
            if f1 > err_u1
                err_u1 = f1;
            end
            if f2 > err_u2
                err_u2 = f2;
            end
            if f3 > err_w3
                err_w3 = f3;
            end
            if f4 > err_T
                err_T = f4;
            end
            if f5 > err_p
                err_p = f5;
            end
        end
        error_u = max(err_u1,err_u2);
        error_w = err_w3;
        error_T = err_T;
        error_p = err_p;
    case "L2"
        err_u1 = 0;
        err_u2 = 0;
        err_w3 = 0;
        err_T = 0;
        err_p = 0;
        for n = 1:mesh.Ne
            V = mesh.node(mesh.elem(n,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "x");
            psi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "x");
            u1h = phi * SOL(fespace1.elem(n,:));
            u2h = phi * SOL(fespace1.elem(n,:)+Nb1);
            w3h = phi * SOL(fespace1.elem(n,:)+2*Nb1);
            Th = phi * SOL(fespace1.elem(n,:)+3*Nb1);
            ph = psi * SOL(fespace2.elem(n,:)+4*Nb1);
            f1 = (pde.u1(int_points_loc(:,1),int_points_loc(:,2),t_end) - u1h).^2;
            f2 = (pde.u2(int_points_loc(:,1),int_points_loc(:,2),t_end) - u2h).^2;
            f3 = (pde.w3(int_points_loc(:,1),int_points_loc(:,2),t_end) - w3h).^2;
            f4 = (pde.T(int_points_loc(:,1),int_points_loc(:,2),t_end) - Th).^2;
            f5 = (pde.p(int_points_loc(:,1),int_points_loc(:,2),t_end) - ph).^2;
            r1 = int_weights_loc * f1;
            r2 = int_weights_loc * f2;
            r3 = int_weights_loc * f3;
            r4 = int_weights_loc * f4;
            r5 = int_weights_loc * f5;
            err_u1 = err_u1 + r1;
            err_u2 = err_u2 + r2;
            err_w3 = err_w3 + r3;
            err_T = err_T + r4;
            err_p = err_p + r5;
        end
        error_u = sqrt(err_u1 + err_u2);
        error_w = sqrt(err_w3);
        error_T = sqrt(err_T);
        error_p = sqrt(err_p);
    case "H1"
        err_u1x = 0;
        err_u1y = 0;
        err_u2x = 0;
        err_u2y = 0;
        err_w3x = 0;
        err_w3y = 0;
        err_Tx = 0;
        err_Ty = 0;
        err_px = 0;
        err_py = 0;
        for n = 1:mesh.Ne
            V = mesh.node(mesh.elem(n,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
            phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "dx");
            phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "dy");
            psix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "dx");
            psiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "dy");
            u1hx = phix * SOL(fespace1.elem(n,:));
            u1hy = phiy * SOL(fespace1.elem(n,:));
            u2hx = phix * SOL(fespace1.elem(n,:)+Nb1);
            u2hy = phiy * SOL(fespace1.elem(n,:)+Nb1);
            w3hx = phix * SOL(fespace1.elem(n,:)+2*Nb1);
            w3hy = phiy * SOL(fespace1.elem(n,:)+2*Nb1);
            Thx = phix * SOL(fespace1.elem(n,:)+3*Nb1);
            Thy = phiy * SOL(fespace1.elem(n,:)+3*Nb1);
            phx = psix * SOL(fespace2.elem(n,:)+4*Nb1);
            phy = psiy * SOL(fespace2.elem(n,:)+4*Nb1);
            f1x = (pde.u1dx(int_points_loc(:,1),int_points_loc(:,2),t_end) - u1hx).^2;
            f1y = (pde.u1dy(int_points_loc(:,1),int_points_loc(:,2),t_end) - u1hy).^2;
            f2x = (pde.u2dx(int_points_loc(:,1),int_points_loc(:,2),t_end) - u2hx).^2;
            f2y = (pde.u2dy(int_points_loc(:,1),int_points_loc(:,2),t_end) - u2hy).^2;
            f3x = (pde.w3dx(int_points_loc(:,1),int_points_loc(:,2),t_end) - w3hx).^2;
            f3y = (pde.w3dy(int_points_loc(:,1),int_points_loc(:,2),t_end) - w3hy).^2;
            f4x = (pde.Tdx(int_points_loc(:,1),int_points_loc(:,2),t_end) - Thx).^2;
            f4y = (pde.Tdy(int_points_loc(:,1),int_points_loc(:,2),t_end) - Thy).^2;
            f5x = (pde.pdx(int_points_loc(:,1),int_points_loc(:,2),t_end) - phx).^2;
            f5y = (pde.pdy(int_points_loc(:,1),int_points_loc(:,2),t_end) - phy).^2;
            r1x = int_weights_loc * f1x;
            r1y = int_weights_loc * f1y;
            r2x = int_weights_loc * f2x;
            r2y = int_weights_loc * f2y;
            r3x = int_weights_loc * f3x;
            r3y = int_weights_loc * f3y;
            r4x = int_weights_loc * f4x;
            r4y = int_weights_loc * f4y;
            r5x = int_weights_loc * f5x;
            r5y = int_weights_loc * f5y;
            err_u1x = err_u1x + r1x;
            err_u1y = err_u2y + r1y;
            err_u2x = err_u2x + r2x;
            err_u2y = err_u2y + r2y;
            err_w3x = err_w3x + r3x;
            err_w3y = err_w3y + r3y;
            err_Tx = err_Tx + r4x;
            err_Ty = err_Ty + r4y;
            err_px = err_px + r5x;
            err_py = err_py + r5y;
        end
        error_u = sqrt(err_u1x + err_u1y + err_u2x + err_u2y);
        error_w = sqrt(err_w3x + err_w3y);
        error_T = sqrt(err_Tx + err_Ty);
        error_p = sqrt(err_px + err_py);
    otherwise
        error("Invalid error type.");
end
end